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INTRODUCTION 


The  purpose  of  this  project  is  to  develop  a  computer  program  that 
will  design  and  implement  a  recursive  digital  filter  to  meet  user  speci¬ 
fications  given  as  inputs  to  the  program.  The  filter  type  is  selected 
from  three  well-known  analog  filter  approximations:  (1)  Butterworth , 

maximally  flat;  (2)  Chebyshev,  equiripple  in  the  passband;  and  (3) 
elliptic  (or  Cauer)  ,  equiripple  in  both  the  passband  and  the  stop 
band.  For  a  low-pass  filter,  for  example,  the  user  specifies  the  filter 
type,  the  filter  order,  the  cutoff  frequency,  and  the  sampling  rate. 
High-pass,  band-pass,  and  band-stop  filters  also  may  be  chosen  for 
design. 

The  program  calculates  the  poles  and  the  zeros  of  the  low-pass 
filter  in  the  analog  s  domain.  The  bilinear  transformation  technique  is 
then  used  to  determine  the  coefficients  of  the  transfer  function  in  the 
z  domain.  If  a  high-pass,  band-pass,  or  band-stop  filter  is  required, 
then  the  coefficients  are  determined  by  a  transformation  in  the  z  domain 
of  the  prototype  low-pass  digital  filter  whose  coefficients  have  been 
determined  from  an  analog  low-pass  filter  using  the  bilinear  transforma¬ 
tion. 


After  determining  the  coefficients  of  the  digital  filter  transfer 
function,  a  plot  of  the  magnitude  and  the  phase  may  be  obtained  if 
desired  by  the  user.  A  list  of  the  coefficients  also  is  printed. 
Another  option  is  to  implement  the  *■  '  ter  as  a  system  of  linear  differ¬ 
ence  equations  and  to  process  input  data  consisting  of  sample  values  at 
any  given  sample  rate  through  the  filter.  The  filtered  output  data 
points  are  then  printed  as  output. 

In  summary,  the  program  performs  the  following: 

a.  Calculates  and  prints  the  coefficients  of  the  desired  digital 
f i 1  ter . 

b.  Plots  the  magnitude  and  phase  responses. 

c.  Implements  the  digital  filter  and  processes  data  through  the 
filter . 

d.  Outputs  the  filtered  data. 


2 .  APPROACH 

There  are  two  general  types  of  digital  filters,  recursive  and  non¬ 
recursive.  The  output  response  of  a  nonrecursive  filter  is  a  function 


of  only  the  present  and  past  values  of  the  input  excitation.  A  recur¬ 
sive  filter  is  one  in  which  the  present  output  response  is  a  function  of 
the  present  and  past  values  of  the  input ,  as  well  as  past  values  of  the 
output . 

As  in  analog  filters,  the  approximation  step  in  the  design  of  digi¬ 
tal  filters  is  the  process  whereby  a  realizable  transfer  function  satis¬ 
fying  prescribed  specifications  is  obtained.  To  be  realizable  as  a 
recursive  filter,  a  transfer  function  must  satisfy  the  following  con¬ 
straints  . 


a.  It  must  be  a  rational  function  of  z  with  real  coeff icients. 

b.  Its  poles  must  lie  within  the  unit  circle  of  the  z  plane. 

c.  The  degree  of  the  numerator  polynomial  must  be  equal  to  or  less 
than  that  of  the  denominator  polynomial. 

Recursive  filter  approximations  can  be  obtained  from  analog  filter 
approximations  by  using  the  following  methods,  all  of  which  satisfy  the 
above  constraints: 

Invariant- impulse-response  method 
Matched  z  transformation 
Bilinear  transformation 

Each  method  has  certain  advantages  and  disadvantages. 

2. 1  Invariant- Impulse-Response  Method 

Given  an  analog  filter  transfer  function,  HA(s),  the  invariant- 
impulse-response  method  is  implemented  as  follows: 

a.  Obtain  the  impulse  response  of  analog  filter  H,(t). 

b.  Replace  t  by  nT  in  Hft(t). 

c.  Form  the  z  transform  of  Hft(nT).  This  gives  HD(z),  the 
digital  filter  transfer  function.  T  is  the  sample  period  =  1/fg. 

This  method  gives  good  results  if  HA(jw)  »  0  for  w  >  u>g/2. 
However,  aliasing  errors  tend  to  restrict  this  method  to  the  design  of 
all  pole  filters. 

2.2  Matched  z  Transformation 


Given  continuous-time  transfer  function 


Ho  n  (s  -  si) 

H  (s)  =  _± - 

A  N 

.11  (S  "  Pi) 

1  =  1 


a  corresponding  digital  transfer  function  can  be  formed  as 


Vz) 


(2  +  1)L 


where  L  is  an  integer  equal  to  the  number  of  zeros  at  s  =  00  in  Hft(s). 
This  method  gives  reasonable  results  for  high-pass  and  band-stop  fil¬ 
ters,  although  it  tends  to  distort  the  passb.tnd  ripple  in  Chebyshev  and 
elliptic  filters.  For  low-pass  and  band-pass  filters,  better  approxima¬ 
tions  can  be  obtained  by  using  the  modified  invariant-impulse-response 
method. 1 


2.3  Bilinear  Transformation 


The  bilinear  transformation  method  yields  a  digital  filter  with 
approximately  the  same  time-domain  response  as  the  original  analog 
filter  for  any  excitation. 


H  ( z)  =  H  (s> 
D  A 


s 


2  z  -  1 
T  z  +  1 


The  bilinear 

unit  circle 

b. 

c . 


transformation  maps  these: 

The  open  right-half  s  plane  onto  the  region  exterior  to  the 
|z|  =  1  of  the  z  plane 

The  j  axis  of  the  s  plane  onto  the  unit  circle  |z|  =1 
The  open  left-half  s  plane  onto  the  interior  of  the  unit 


circle 


From  property  c  it  follows  that  a  stable  analog  filter  yields  a 
stable  digital  filter  and,  since  the  transformation  has  real  coeffi¬ 
cients,  Hd(z)  has  real  coefficients. 


lA.  Antoniou,  Digital  Filters:  Analysis  and  Design,  McGraw-Hill  Book 
Co.,  Inc.,  New  York  (1979). 
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If  only  the  amplitude  response  is  of  concern,  the  warping 
effect  can  for  all  practical  purposes  be  eliminated  by  prewarping  the 
analog  filter.  For  example,  if  a  low-pass  cutoff  frequency  of  is 

desired  for  the  digital  filter,  then  the  analog  filter  is  first  designed 
for  an  unnormalized  cutoff  frequency  given  by 


The  phase  response  of  the  derived  digital  filter  is  nonlinear 
because  of  the  warping  effect.  Furthermore ,  little  can  be  done  to 
linearize  it  except  by  employing  delay  equalization.  Consequently,  if 
it  is  mandatory  to  preserve  a  linear  phase  response,  alternative  methods 
should  be  considered. 


The  bilinear  transformation  is  the  most  important  of  the  tech¬ 
niques  used  to  obtain  digital  recursive  filters  from  analog  filters. 
The  passband  and  loss  character  istics  of  the  analog  filter  are  pre¬ 
served,  and  there  is  no  aliasing  effect.  Frequency  distortion  is  com¬ 
pensated  at  and  the  transfer  function  can  be  obtained  by  a 

relatively  easy  transformation.  For  these  reasons,  the  bilinear  trans¬ 
formation  was  chosen  as  the  method  of  obtaining  the  digital  filter 
transfer  function  for  this  project. 


2.4  Analog  Filter  Designs 


The  theory  of  analog  filter  approximations  has  been  extensively 
developed. ^  4 


2. 4. 1  Butterworth 


The  Butterworth  approximation  is  the  simplest  all-pole  type. 
In  the  stop  band,  the  attenuation  characteristics  monotonically  decrease 
as  a  function  of  frequency.  The  Butterworth  filter  is  termed  also  a 
maximally  flat  magnitude  approximation  in  that  the  error  in  the  passband 
is  also  a  monotonically  decreasing  function.  The  equations  for  calcula¬ 
ting  the  pole  locations  are  given  in  appendix  A. 


Antoniou,  Digital  Filters:  Analysis  and  Design,  McGraw-Hill  Book 
Co.,  Inc.,  New  York  (1979). 

2M.  S.  Ghausi,  Principles  and  Design  of  Linear  Active  Circuits, 
McGraw-Hill  Book  Co.,  Inc.,  New  York  (196 5),  ch.  4. 

2H.  Y.-F.  Lam,  Analog  and  Digital  Filters,  Prentice-Hall,  Inc.,  Engle¬ 
wood  Cliffs,  NJ  (1979). 

4 D .  E.  Johnson,  Introduction  to  Filter  Theory,  Prentice-Hall ,  Inc., 
Englewood  Cliffs,  NJ  (1976). 
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2.4.2  Chebyshev 

A  second  approximation  that  improves  on  the  rate  of  change  of 
the  attenuation  between  passband  and  stop  band  over  that  of  the  Butter- 
worth  filter  is  the  Chebyshev  filter.  The  error  in  the  passband  is 
distributed  evenly  in  an  oscillating  manner.  This  is  called  an  equi- 
ripple  approximation.  In  the  stop  band,  the  magnitude  decreases  mono- 
tonically  with  a  faster  cutoff  rate  than  that  of  a  Butterworth  filter  of 
the  same  order.  The  Chebyshev  filter  is  an  all-pole  type  and  is  based 
on  the  nth-order  Chebyshev  polynomials,  Cn(u>)  • 


cos  (n  cos”1  uj)  ,  0  <  u>  <  1  , 
cosh  (n  cosh”1  j)  ,  u>  >  1  . 


The  recursion  formula  for  finding  the  nth-order  polynomial  is 

CQ  (oi )  =  1  , 

C 1  (cj )  =io  , 

C  (u>)  =  2u>C  .(to)  -  C  „{w)  . 
n  n-1  n-2 

The  equations  for  calculating  the  pole  locations  are  given  in  appen¬ 
dix  A. 


2.4.3  Elliptic 

The  third  type  of  filter  approximation  considered  in  this 
report  is  called  the  elliptic  approximation  and  was  first  introduced  by 
Wilhelm  Cauer.  This  approximation  is  based  on  the  Jfecobi  elliptic  sine 
functions.  In  this  approximation,  the  error  in  the  passband  is  again 
distributed  evenly  in  an  oscillating  manner.  However,  instead  of  a 
monotonically  decreasing  characteristic  in  the  stop  band,  the  stop-band 
attenuation  oscillates  between  infinity  and  a  prescribed  maximum.  Thus, 
there  is  an  equiripple  characteristic  in  both  passband  and  stop  band. 
The  elliptic  approximation  is  more  efficient  than  the  Butterworth  and 
Chebyshev  approximations  in  that  the  transition  between  passband  and 
stop  band  is  steeper  for  a  given  filter  order.  An  elliptic  filter 
transfer  function  has  both  poles  and  zeros  and  has  been  shown  to  be 
optimal  in  the  sense  of  having  the  sharpest  transition  of  any  approxi¬ 
mation.  The  technique  used  to  calculate  the  coefficients  for  the  ellip¬ 
tic  filter  is  given  in  appendix  A  and  was  taken  from  the  development  by 
Antoniou.1 


1A.  Antoniou,  Digital  Filters:  Analysis  and  Design,  McGraw-Hill  Book 
Co.,  Inc.,  New  York  (1979). 
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2.5  Digital  Filter  Realization 


There  are  two  possible  canonical  forms  of  linear  difference 
equations  that  can  be  realized  from  the  z  domain  transfer  function, 
H(z).  One,  the  cascade  form,  follows  from  the  factored  form  of  H(z). 
The  other,  the  parallel  form,  requires  the  expansion  of  H(z)  into 
partial  fractions.  The  parallel  form  was  chosen  for  this  project 
because  of  its  reduced  sensitivity  to  noise. 

Realization  of  a  recursive  digital  filter  in  parallel  form  is 
shown  in  figure  2. 


INPUT 


OUTPUT 


Z-V,.  +  Ait.  z'1 

- 2! - li - 

k=  1  1  +  ®1k  r'  +  ®2k  r* 


Figure  2.  Recursive  digital  filter  in  parallel  form. 


2.6  Band  Transformations 


For  this  project,  high-pass,  band-pass,  and  band-stop  filters 
are  obtained  from  the  digital  low-pass  prototype  by  a  transformation  in 
the  z  domain.  Another  approach  is  to  design  the  analog  low-pass  proto¬ 
type,  perform  the  transformation  in  the  s  domain  to  high  pass,  band 
pass,  or  band  stop,  and  then  transform  the  resulting  transfer  function 
to  the  z  domain  by  the  bilinear  transformation. 

The  work  of  Constantinides5  was  used  to  develop  the  z  domain 
transformations.  The  equations  are  given  in  appendix  A. 


5A.  G.  Constantinides,  Spectral  Transformations  for  Digital  Filters, 
Proc.  IEE,  U7_  ( August  1970),  1585-1590. 


3.  FORTRAN  IV  PROGRAM 

3. 1  Outl ine 

A  FORTRAN  IV  computer  program  was  written  using  the  equations 
and  the  techniques  discussed  in  section  2.  The  program  was  designed  in 
a  modular  form  to  allow  for  easy  modification  as  required.  Hie  flow 
chart  of  the  program  is  given  in  figure  3.  The  main  part  of  the  program 
is  simply  a  controller  that  calls  in  subroutines  as  required.  The 
blocks  in  a  column  in  the  center  of  the  flow  chart  are  all  in  the  main 
program,  and  the  blocks  to  the  right  or  the  left  of  center  are  separate 
subroutines.  For  instance,  if  a  different  form  of  input  data  is  re¬ 
quired  or  a  different  output  format  is  desired,  then  only  the  subroutine 
dealing  specifically  with  that  function  need  be  modified.  If  a  differ¬ 
ent  analog  design  type  or  technique  is  to  be  used,  it  can  be  either 
substituted  for  one  of  the  existing  three  design  types  or  added  on  by 
extending  a  "computed  go  to”  statement  in  the  main  program. 

All  data  are  stored  in  labeled  and  unlabeled  common  areas  of  a 
subroutine  simply  by  including  the  proper  common  statement. 

This  program  was  developed  by  using  an  IBM  System/370  com¬ 
puter.  The  H  Extended  Plus  FORTRAN  compiler  was  used  with  the  auto¬ 
double  option  specified.  This  option  doubles  the  overall  precision 
specified  in  the  program.  Maximum  possible  precision  is  used  for  all 
the  design  calculations  as  well  as  in  the  implementation  itself.  This 
precision  is  16  bytes  for  real  variables  and  32  bytes  for  complex  varia¬ 
bles.  The  complete  program  listing  is  given  in  appendix  B. 

The  flow  of  the  program  is  straightforward  because  it  begins  at 
the  top  and  progresses  without  any  diversions  to  the  bottom  of  the 
chart.  The  basic  steps  of  the  program  and  the  general  function  of  each 
subroutine  are  as  follows. 

3.1.1  RDATA 

After  the  program  is  initiated,  it  immediately  calls  the 
RDATA  subroutine.  This  subroutine  reads  in  the  necessary  filter  design 
criteria  and  the  data  samples  to  be  processed.  Only  batch  processing  of 
data  is  possible  with  the  existing  program,  but  the  conversion  to  real 
time  processing  would  not  be  difficult.  The  input  data  array  is  dimen¬ 
sioned  for  up  to  1000  amplitude  points  at  the  specified  sampling  fre¬ 
quency.  If  only  a  filter  design  is  required,  then  the  number  of  data 
points  (NTIMES)  is  set  equal  to  zero,  and  RDATA  bypasses  the  input  of 
data  to  be  processed.  Setting  NTIMES  equal  to  zero  bypasses  also  the 
implementation  subroutine.  RDATA  also  prints  out  the  design  criteria 
data  for  future  reference. 

The  specific  format  required  for  data  input  is  given  in 
section  3.2. 
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3.1.2  WARP 


The  RDATA  subroutine  returns  control  to  the  main  proqram, 
which  calculates  the  prewarped  analog  design  frequency  (WARP)  and  then 
calls  the  appropriate  analog  design  subroutine.  The  analog  filter 
subroutines  calculate  the  poles  and  the  zeros  (  if  any)  of  the  squared 
magnitude  transfer  function  that  '  ie  in  the  left  half  of  the  complex  s 
plane.  The  resulting  transfer  function  aff-"»r  the  right-half  plane  poles 
and  zeros  are  eliminated  is  of  the  f*  im 

( AMP)  (s  -  z^s  )  .  .  .  U  -  z^) 

H(  s)  - - - , 

(S  “  Pl)(S  -  P2)  *  ’  (S  “  Pi) 

where  z^  and  are  complex  conjugates  of  z^  and  p^  (for  all  i).  The 
equations  needed  to  calculate  these  poles  and  zeros  are  given  in  appen¬ 
dix  A  and  are  derived  as  discussed  in  section  2. 

3.1.3  FACTOR 


The  FACTOR  subroutine  is  then  called.  FACTOR  expands  the 
transfer  function  into  partial  fractions.  If  the  numbers  of  poles  and 
zeros  are  equal,  then  this  expansion  has  a  constant  gain  equal  to  the 
amplitude  multiplier  (AMP)  of  the  transfer  function.  FACTOR  then  com¬ 
bines  the  complex  conjugate  pairs  and  applies  the  bilinear  transfor¬ 
mation  , 


to  obtain  the  coefficients  for  the  parallel  second-order  sections  of  the 
digital  filter.  The  general  form  of  the  resulting  z  domain  transfer 
function  is^ 


H(z~l) 


(--->)  f , >oi  *  _2 , 

(—>  1  +  BUZ_1  +  B2pZ  z 

i=1 


where  r  =  n/2  for  n  even,  r  =  (n  +  1)/2  for  n  odd,  and  n  is  the  order  of 
the  filter.  The  coefficients  A,^,  A^,  B^,  and  B2^  are  calculated  by 
FACTOR  using  the  analog  poles  ana  zeros.  If  the  order  of  the  filter  is 
odd,  then  there  is  one  real  pole,  so  the  A^  and  B2  coefficients  are 
equal  to  zero  for  that  pole. 

^G.  C.  Temes  and  S.  K.  Mitra ,  Modern  Filter  Theory  and  Design,  John 
Wiley  S  Sons,  Inc.,  New  York  (1973). 
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3.1.4 


F  PLOT 


A  plot  of  the  frequency  resj^onse  (maqnituiie  and  phase)  of 
the  low-pass  diqital  filter  is  then  made  by  the  FPLOT  subroutine.  It  is 
plotted  by  letting 

_i  -jiitT  m  m 

z  ‘  =  o  -  cos  i.iT  -  •)  sin  oT 

in  the  z  domain  transfer  function  given  above. 

The  plots  given  in  this  report  were  made  on  a  Houston 
Instruments  plotter  using  existing  plotting  software.'’  Both  the 
amplitude  and  phase  responses  are  plotted  over  the  normalized  range  of 
frequencies  of  zero  to  one.  That  is,  the  frequency  scale  is  normalized 
to  the  sampling  frequency.  FPLOT  also  prints  the  values  of  the  digital 
filter  coefficients  to  25  decimal  places  for  •‘uture  use. 

3.1.5  Filter  Choice 


The  main  proqram  then  determines  if  a  high-pass,  band-pass, 
or  band-stop  filter  is  desired  or  if  the  low-pass  filter  is  ready  to  be 
implemented . 

The  transfer  function  for  the  required  high-pass,  band-pass, 
or  band-stop  filter  is  obtained  by  r  appropriate  trans f ormat ion  in  the 
z  domain  of  the  low-pass  digital  filter.  Band-pass  and  band-stop  trans¬ 
formations  result  in  a  doublinq  of  the  filter  order. 

After  tr  u.sformation ,  the  frequency  characteristics  of  the 
new  filter  are  plotted  by  FPLOT. 

3 .  1 .  f .  FILIMP 


The  actual  implementation  of  the  digital  filter  is  simple. 
The  FILIMP  subroutine  performs  this  function  and  processes  the  data  to 
be  filtered.  Application  of  the  inverse  z  transform  to  the  general  form 
for  the  z  domain  transfer  function  results  in  the  recursive  difference 
equation 

Yout(kT)  =  X(kT)  +  X[(k  -  1)T1  , 


' Thomas  V.  Noon  and  Egon  Marx,  User's  Manual  for  the  Modular  Analysis- 
Package  Libraries  ANAPAC  and  TRANL,  Harry  Diamond  Laboratories  HDL-TR- 
1782 -S  (September  1978). 
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where 


r 

X(kT)  =  l  Xi(kT)  , 
i=1 

Xi(kT)  =  AQiYin(kT)  +  AuYin[(k  -  1  ) T ] 

-  BuXi[(k  -  1)T]  -  B2iXi[(k  -  2)T]  • 

3.1.7  PDATA 

The  PDATA  subroutine  provides  plots  of  the  input  and  output 
data.  This  routine  may  be  easily  modified  to  provide  the  output  data  in 
any  desired  format. 

3.2  Program  Use 

There  are  at  most  nine  parameters  that  must  be  supplied  to  the 
program  for  the  design  of  a  filter.  The  minimum  number  needed  is  five 
for  the  design  of  a  Butterworth  low-pass  filter.  This  minimum  set, 
which  is  used  in  all  designs,  consists  of  these  parameters: 

a.  1TYPE,  the  type  of  analog  design 

1  =  Butterworth 

2  =  Chebyshev 

3  =  Elliptic  (Cauer) 

b.  N,  the  order  of  the  filter 

c.  ITRANS,  corresponding  to  the  band  transformation  desired 
0  =  None 

1  =  High-pass  transformation 

2  =  Band-pass  transformation 

3  =  Band-stop  transformation 

d.  FC,  the  cutoff  frequency  desired  ( in  hertz) 

e.  FS,  the  sampling  frequency  of  the  data  (in  hertz)  to 

which  the  filter  will  be  normalized 
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If  a  low-pass  Chebyshev  is  required,  then  an  additional  param¬ 
eter,  EPS  1 ,  the  minimum  allowed  amplitude  in  the  passband,  also  must  be 
entered.  Similarly,  for  elliptic  filters,  EPS  1  must  be  given  along  with 
EPS2,  the  transition  region  selectivity  factor.  The  maximum  value  for 
EPS2  is  0.95.  All  of  these  parameters  are  read  from  the  same  data 
card.  These  are  the  formats  for  all  the  data  cards  and  the  order  in 
which  the  cards  are  read  by  the  program : 

ITY PE , N , ITRANS , FC , FS , E  PS  1 , EPS  2 
Format — 315,  4E10.3 

If  ITRANS  =  0,  skip  reading  FI  and  F2. 

F 1 ,  F2  , 

Format — 2E10. 3 

NTIMES 
Format — 15 

If  ITIMES  =  0,  skip  reading  VIN. 

VIN 

Format — 8E10.3 — maximum  of  125  cards 
Heading  for  low-pass  design 
If  ITRANS  =0,  skip  next  heading. 

Heading  for  band  transformation 
If  NTIMES  =  0,  skip  next  two  headings. 

Heading  for  VIN  plot 
Heading  for  VOUT  plot 

If  a  band  transformation  is  specified,  then  an  additional  data 
card  is  needed  that  has  the  lower  cutoff  frequency,  FI  (also  the  high- 
pass  cutoff  frequency) ,  and  the  upper  cutoff  frequency,  F2,  for  the 
band-pass  and  band-stop  filters. 

Next,  the  number  of  data  sample  points  is  read.  If  NTIMES  is  set 
equal  to  zero,  then  no  input  data  are  read,  and  the  program  is  termi¬ 
nated  after  completion  of  the  filter  design. 

Separate  titles  are  read  in  for  the  low-pass  filter  transfer 
function  plots,  any  band  transformation  frequency  plots,  and  the  heading 
to  appear  on  the  data  plots.  These  titles  can  be  up  to  80  characters  of 
any  form  desired. 
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A  sample  filter  design  is  given  now  to  illustrate  the  order  and 
the  manner  in  which  the  data  are  assembled.  A  fourth-order  elliptic 
band-pass  filter  is  chosen  since  this  design  requires  input  of  all  the 
design  parameters.  We  have  this: 

ITYPE  =  3 

N  =4 

ITRANS  =  2 

Now  let  us  assume  that  the  sampling  frequency  is  10,000  Hz  and  that  the 
passband  is  to  be  from  1000  to  3000  Hz.  In  addition,  we  allow  the 
passband  amplitude  to  vary  between  1.0  and  0.9  and  minimize  the  transi¬ 
tion  between  the  passband  and  the  stop  band.  This  action  means  that  the 
cutoff  frequency  for  the  low-pass  design  must  be  2000  Hz,  which  is  the 
width  of  the  passband.  Therefore,  the  remaining  parameters  are  these: 

FC  =  2000. 0 

FS  =  10000.0 

EP SI  =0.9 

EPS2  =  0.95 

FI  =  1000.0 

F2  =  3000.0 

These  data  are  given  in  figure  4  as  they  appear  on  the  input  data 
cards.  The  plot  title  cards  also  are  shown  in  the  figure.  If  only  a 
low-pass  design  is  requested,  then  the  second  title  card  is  omitted. 
Similarly,  if  no  data  points  are  given,  then  the  third  title  card  as 
well  as  the  input  data  cards  are  omitted. 

The  resulting  printout  for  this  example  problem  is  given  in 
figure  5.  The  low-pass  transfer  function  is  given  in  figures  6  and  7, 
and  the  band-pass  function  is  given  in  figures  8  and  9.  Running  this 
example  required  6.85  s  in  the  central  processing  unit  (CPU)  to  compile 
the  FORTRAN  coding,  2.44  s  in  the  CPU  to  link  and  edit  the  object  mod¬ 
ules  with  library  functions,  and  about  2  s  in  the  CPU  to  perform  the 
actual  calculations.  The  filtering  of  101  data  points  required  less 
than  2  s  of  CPU  time.  Additional  examples  of  filter  designs  are  given 
in  section  4. 
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Figure  4.  Input  data  for  example  1 . 
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Figure  5.  Output  data  for  example  1. 


Figure  6.  Low-pass  Cauer  filter:  amplitude  versus  frequency 


Figure  7.  Low-pass  Cauer  filter:  phase  versus  frequency 
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4.  DESIGN  EXAMPLES 


Several  examples  are  given  in  this  section  to  illustrate  the  differ¬ 
ences  in  the  various  analog  designs  and  also  to  show  the  results  of  band 
transformations.  In  addition,  a  simple  waveform  consisting  of  three 
sinusoids  at  different  frequencies  is  used  to  demonstrate  the  filtering 
action  of  the  various  digital  filters.  . 

4. 1  Analog  Design 

The  three  types  of  analog  designs  are  discussed  in  section  2. 
A  graphical  presentation  of  the  basic  differences  in  these  three  designs 
is  given  here.  A  fourth-order  elliptic  low-pass  design  is  given  in 
section  3.2.  Similar  fourth-order  digital  filters  using  the  Butterworth 
and  Chebyshev  analog  designs  are  presented  in  figures  10  to  13.  Compar¬ 
ing  these  transfer  functions,  we  can  easily  see  that  the  Butterworth 
filter  offers  the  smoothest  response  in  both  the  passband  and  the  stop 
band.  The  sharper  transition  from  passband  to  stop  band  of  the 
Chebyshev  design  over  the  Butterworth  is  obtained  at  the  expense  of 
smoothness  in  the  passband.  The  elliptic  design  obviously  has  the 
sharpest  transition  of  the  three  designs.  However,  the  ripple  in  both 
the  passband  and  the  stop  band  is  the  penalty  for  this  sharp  transi¬ 
tion.  The  design  to  use  is  chosen  according  to  the  requirements  of  the 
specific  filtering  problem  being  considered. 

4. 2  Band  Transformations 

An  example  of  a  band-pass  transformation  is  given  in  figures  8 
and  9.  The  low-pass  Butterworth  filter  of  figures  10  and  11  was  trans¬ 
formed  to  a  high-pass  filter  with  a  cutoff  of  3000  Hz  as  shown  in  fig¬ 
ures  14  and  15.  The  band-stop  dual  of  the  passband  filter  was  made  by 
using  the  Chebyshev  design  and  is  shown  in  figures  16  and  17. 
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Figure  14.  High-pass  Butterworth  filter:  amplitude  versus  frequency 
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Figure  15.  High-pass  Butterworth  filter:  phase  versus  frequency 
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Figure  16.  Band-stop  Chebyshev  filter:  amplitude  versus  frequency 


4. 3  Data  Processing 


A  simple  waveform  was  generated  to  illustrate  the  operation  of 
the  filters  on  input  data.  The  waveform  used  was  the  sum  of  three 

sinusoids  at  0.05,  0.15,  and  0.35  times  the  sampling  frequency  with  peak 
amplitudes  of  3,  2,  and  5,  respectively.  This  input  waveform  is  shown 

in  figure  18.  A  20th-order  Butterworth  low-pass  design  was  used  to 
eliminate  the  two  highest  frequency  sine  waves,  as  shown  in  figure  19. 
A  filter  of  high  order  was  used  to  demonstrate  the  group  delay  effect  of 
the  filter.  This  effect  shows  at  the  beginning  of  the  output  waveform 
as  a  delay  time  before  any  output  occurs.  Next,  a  fourth-order  Butter- 
worth  high-pass  filter  was  used  to  eliminate  the  two  low-frequency 

signals.  This  result  is  shown  in  figure  20.  The  low  sampling  rate 

caused  the  distortion  on  the  sine  wave  seen  in  this  figure.  The  initial 
design  example  was  used  to  pass  only  the  center  sinusoid  as  shown  in 
figure  21.  Here,  again,  the  effect  of  the  low  sampling  rate  is  seen  as 
a  distortion  of  the  waveform.  However,  the  effect  is  not  as  great 

because  of  the  lower  frequency  of  the  sinusoid.  The  Chebyshev  band-stop 
dual  of  this  passband  filter  was  then  used  to  eliminate  the  middle  sine 
wave.  This  waveform  is  given  in  figure  22. 
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Figure  18.  Input  waveform  obtained  by  summing  three  sinusoids. 
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Figure  19.  Output  waveform  obtained  by  using  low-pass  Butterworth 
filter . 


Figure  20 . 


Output  waveform  obtained  by  using  high-pass  Butterworth 
filter. 
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Figure  21.  Output  waveform  obtained  by  using  band-pass  Cauer  filter 
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Figure  22.  Output  waveform  obtained  by  using  band-stop  Chebyshev 
filter. 


5.  CONCLUSIONS  AND  RECOMMENDATIONS 


A  versatile  program  has  been  developed  for  designing  and  implemen¬ 
ting  recursive  digital  filters.  Three  well-known  analog  filter  designs 
were  used  with  the  bilinear  transformation  to  obtain  the  desired  digital 
filters.  The  bilinear  technique  was  chosen  because  it  eliminates  the 
effects  of  aliasing  that  occur  with  other  techniques.  However,  the 
resulting  warping  of  the  frequency  scale  may  not  be  acceptable  in  some 
applications.  The  extent  to  which  this  warping  affects  the  resulting 
output  was  not  investigated.  Examples  using  the  different  analog  de¬ 
signs  are  given  in  section  4  along  with  transformations  from  low-pass  to 
high-pass,  band-pass,  and  band-stop  designs.  All  of  the  filter  design 
calculations  as  well  as  the  implementations  use  the  highest  precision 
possible  for  the  IBM  System/370  computer.  No  investigation  was  made  to 
determine  the  effects  of  using  lower  precision. 

All  of  the  equations  used  in  developing  the  code  are  included  in 
appendix  A,  and  the  complete  computer  program  code  is  listed  in  appendix 
B.  A  bibliography  of  related  publications  also  is  included.  This 
bibliography  is  by  no  means  complete  since  the  amount  of  information 
published  seems  almost  endless.  Although  there  is  a  large  amount  of 
information  available,  no  one  source  proved  adequate  in  presenting  all 
of  the  steps  necessary  to  completely  develop  a  digital  filter.  It  is 
hoped,  therefore,  that  this  document  will  prove  helpful  to  others  ven¬ 
turing  into  this  field  for  the  first  time. 

Although  the  program  developed  here  is  completely  operable,  there 
are  many  ways  in  which  the  program  may  be  improved  or  modified  to  better 
suit  a  particular  filtering  need.  Also,  investigation  and  analysis  of 
the  various  errors  associated  with  the  techniques  used  should  be  con¬ 
ducted  . 
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APPENDIX  A 


The  following  equations  were  used  in  the  development  of  the  com¬ 
puter  code  for  designing  and  implementing  digital  filters.  All  of  the 
equations  are  referenced  in  the  main  body  of  this  report.  The  applica¬ 
tion  of  each  set  of  equations  is  given  in  the  flow  chart  of  the  program 
(fig.  3). 

A-1.  PREWARPING  (WARP) 


A- 2.  8UTTERW0RTH 

The  transfer  function  for  a  Butterworth  low-pass  filter  of  order  n 
is  given  by 


where 


H_ 


H(s)  = 


0 


",  (■  - 


H„  =  ui  , 
0  c 


=  cutoff  frequency, 

~  j[(*/2)+u(2k-1)/2n]  ,  ,  . 

Pk  =  uce  '  k  =  2' 


n . 


A- 3.  CHEBYSHEV 

The  squared  magnitude  function  is  of  the  form 


|H( ju) |2 


1 

1  +  e2C2(w) 


where 


e  =  ripple  factor, 

Cn(w)  =  Chebyshev  polynomial. 
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For  the  program,  e  is  entered  as 


EPS1  =  1 


(1  +  e2)]/2 

The  transfer  function  is  of  the  form 


H(s)  = 


where 


n  (s  -  Pk) 

k=1 


Pk  =  sin  sinh  (v)  +  j  cos  cosh  (v)  = 

(2k  -  1),  k  =  n  +  1,  n  +  2,  .  .  2n , 


v  =  l  sinh"'  (-), 


(n-1 )/2 


ta(  n+ 1 ) /2~)  JJ,  (ak  +  pk)'  n  odd' 


H0  = 


n/2 


( 177  ,n  H  +  &2y)'  “  «"“■ 


k=  t 
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ELLIPTIC  (CAUER) 


For  the  normalized  filter  u)c  =  1, 
n  =  filter  order, 


k  =  selectivity  factor  =  ujp/u)a, 

Ap  =  maximum  passband  attenuation  (in  decibels), 
Ag  =  minimum  stop-band  attenuation  (in  decibels). 


Of  the  parameters  Ap,  Afl,  k,  and  n,  only  three  are  independent. 
If  three  of  the  four  are  specified,  the  fourth  is  automatically  fixed. 
For  this  project,  n,  k,  and  Ap  are  specified: 

k  =  EPS  2, 

Ap  =  20  log  ( 1  -  EPS1) . 


Aa  is  then  given  by 


A  =10  log 

3 


-0.  1 A 

10  P  -  1 


+  1 


16qr 
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H(s) 


\ 


\ ( n+1 ) /2  ,  n  odd  , 
I n/2 ,  n  even  , 


s  =  2  1 
s  T 


-  z  ' 


1  +  z 


L 

H(z~l)  =  (l  +  z~l)  ^ 


'Ok 


+  A1kz 


-1 


1  +■  B,.z'‘  +  B„,  z 


k=1 


Ik 


2k 


pk = 

_Jk  + 

k  ' 

A0k  = 

M'  * 

V2)  -  MV2)]'0 

A1k  = 

-  DvO 

-  V2)  *  UV2)] 

B1k 

-2[’  -( 

V2)2  -  (V2)2]'0 

B2k  = 

[0  -  \ 

nf  *  (v2)2]/°  • 

D  -  0  +  v2)2  +  (\/2Y  ■ 


A-6.  MAGNITUDE  AND  PHASE  ( F PLOT ) 


The  frequency  response  of  the  designed  digital  filter  is  calcula¬ 
ted  by  substituting  e  -'uiT  for  z  1  in  the  transfer  function  and  then 
using  complex  arithmetic  to  obtain  the  magnitude  and  the  phase  as  a 
function  of  frequency. 


E_1  -  e"3wT  = 


W  =  cos  wT  -  j  sin  wT  , 
r 


k=1 


V  A0k  +  A1kW 
H(W)  =  (1  +  W)  )  -  . 

^  1  *  S1kW  +  B2k«2 
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HIGH-PASS  TRANSFORMATION 


APPENDIX  A 


A  low-pass  digital  filter  of  cutoff  frequency  0  is  first  designed 
for  the  desired  order,  n,  and  type.  Let  ojp  =  the  desired  high-pass 
cutoff  frequency. 

Then  substituting  for  z  ^  in  the  low-pass  transfer  function, 

,  z  *  +  a 


1  +  uz‘ 


where 


(K±) 

(^) 


results  in  the  transfer  function  of  the  desired  high-pass  filter. 


Z  Ok 

771ZTT 


A0k  +  A1kz  * 


+  B2kz  " 


A0k  =  n  -  a) 


A0k  aA  Ik 


A1k  =  (1  "  a) 


aA0k  "  A1k 


B1k  =  [2«  -  B1k(1  +  a2' 


+  2B2kaJ/D'  , 


=  (a2  - 


aB1k  +  B2k 


D'  =  1  -  aB1k  +  a2B2k  . 


A-8.  BAND-PASS  TRANSFORMATION 


A  low-pass  digital  filter  of  cutoff  frequency  6  is  first  designed 
for  the  desired  order  and  type.  The  cutoff  frequency  equals  the  band¬ 
width  of  the  desired  band-pass  filter,  F2  -  FI,  where  F2  and  FI  are  the 
upper  and  lower  cutoff  frequencies  (in  hertz),  respectively. 
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A  substitution  for  z-^  in  the  low-pass  transfer  function  results 
in  the  transfer  function  of  the  desired  band-pass  filter. 

In  general,  this  transformation  is1'2 


>-l 


2uk 
k  +  1 


k  -  1 
k  +  1 


k  +  1 


2otk 
k  +  1 


r-1 


+  1 


where 

cosn  +  F2)T 
u  _  cosir  (F2  -  FAT' 

k  =  cot  [ it ( F2  -  Fi)t]  tan  (it0T). 

By  letting  k  =  1,  then  0  =  F2  -  FI,  and  the  transformation  is  simpli¬ 
fied. 


z 


-L 


The  resulting  band-pass  transfer  function  has  the  form 


n/2 


A11kz  1  +  A01k 


C1  /  11k 


B2 1kz 


-2 


A12kz  1  “  A02k 
1  —  B ^ 2^ z  ^  ^  ^ 


22k 


*A.  G.  Constant inides ,  Spectral  Transformations  for  Digital  Filters, 
Proc.  IEE,  117  (August  1970),  1585-1590. 

2A.  Antoniou,  Digital  Filters:  Analysis  and  Design,  McGraw-Hill  Book 
Co.,  Inc.,  New  York  (1979). 
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for  n  even.  For  n  odd,  there  is  one  real  pole,  and  the  transfer  func¬ 
tion  is  of  the  form 


Aq[1  -  z 


-2  ' 


1  -  a[  1  -  BtJ z"1  -  B ^ : 


-2 


(n-1 )/2 


/  A11kz  *  +  A01k 
^  y1  +  B1 1kz  1  +  B2  Ik2 


(1 


A1 2kz  1  "  A02k 
1  “  B12kz  1  +  B22kz 


A-9.  BAND-STOP  TRANSFORMATION 

The  procedure  for  the  band-stop  transformation  is  the  same  as  for 
the  band-pass  transformation.  In  general,  the  transformation  is 


where 


,-l 


2a  1  -  k 

_2  -  Z  1  +  - 

1  +  k  1  +  k 


1  -  k  2a 

-  z  ^  -  - 

1  +  k  1  +  k 


,-l 


+  1 


COSH^F1  +  F  \T 

a  =  - - - , 

cosw(F2  -  F^T 

k  =  tan  [h(f2  -  Fj)t]  tan  (u0T). 


By  letting  k  =  t,  then  0  =  Fg/2  -  (Fj  -  F.j),  and  the  transforma¬ 
tion  is  simplified. 
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The  resulting  band-stop  transfer  function  has  the  form 


hbs(z_1)  =  l1  -  2az 


-1 


+  z 


-2 


^  /  A11kz  1  "  A01k 

^  1^1  v  "  B'"kZ  1  + 


2  Ik 


,-2 


A12kz  1  +  A02k 
1  +  B12kz  1  +  B22kz  ^ 


for  n  even.  For  n  odd,  the  transfer  function  has  the  form 


hbs(z_1) 


( 1  -  2az-l  +  z“2)Aq 
1  -  a(  1  +  Bjz-1  +  B.,z_2 


+  (  1  -  2az  1  +  z” 


(n-1  )/2  /  -l  .  j,' 

Zl  A11kz  A0 Ik 

K=)  V  -  B11kz_1  +  B2 Ik2 


A1 2kz  1  +  A02k 
1  +  B1 2k z  1  +  B22kz  2' 


A- 10.  FILTER  IMPLEMENTATION 

From  the  z  domain  general  transfer  equation 


V-1  A0i  +  A1iz  1 

H(z'l)  =  (1  +  z-1)  )  - - - - 

Z_ I  l+B-.z1  +  B,  •  z  L- 

i=1  11 


a  time-domain  difference  equation  is  obtained.  The  difference  equation 
is  then  used  to  process  data  samples. 

Yout(kT)  =  X(kT)  +  xHk  -  1)T]  , 
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where 


X(kT) 


Xi(kT) 


l  X^kT), 
i=  1 


AY.  (kT)  +  A.. Y,  [(k  -  1)T] 
Ol  xn  lx  in 


-BuXi[(k  -  1  )T]  -  B2iXi[(k  -  2)T], 
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SUBROUTINE  BAN DP 
SUBROUTINE  BANDS 
SUBROUTINE  FILIMP 
SUBROUTINE  PLT DAT 
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The  main  program  and  all  of  the  subroutines  for  recursive  digital 
filters  are  included  in  this  appendix.  The  only  exceptions  are  standard 
library  functions  and  the  plotting  routine.  If  a  machine  other  than  an 
IBM  is  used,  it  may  be  necessary  to  modify  the  format  and  data  state¬ 
ments  . 


B— 1 •  MAIN  PROGRAM 


REALMS  40*20),  41*20),  B  1*201,  B2(20)»  PI ,A2 
REAl«8  AMP 

CQMPLE  X* 16  P  *2D ) ,  l *20) 

COMMON  FC,  FS,  EPSl,  FPS2,  PI,  1TYPE,  N ,  ) TRANS 
C0MM0N/VDATA/VlNtl000 I,  V0UT  flOOO)  ,NT IMES 
C0MM0N/TRANS/  F 1 , F ? 

COMM ON /F ILT /WC«  W1,W2,AMP,A0,A1,A2,  Bl,  B2,  P,  2 
C  ITVPE-  FILTER-  BUTT ER WORTH 1 1 1 ,  CHEB YSHE VI2 ) ,  ELLIPT1C*3) 

C  N  -  ORDER  OF  FILTER,  20  MAX. 

C  1  TRANS-  NON  E  *  C)  ,  HPI1J,  BP*2),  BSI3) 

C  FC-  COTOFF  FREQUENCY 

C  FS-  SAMPLE  FREQUENCY 

C  EPS1  -  MINIMUM  AMPLITUDE  IN  PASSBAND 

C  EPS2  -  TRANSITION  COEFFICIENT  FOR  CAUER,  0.95  MAX. 

C  FI  -  LOWER  CJT0FF  FREQUENCY,  PASS  AND  STOP  FILTERS  ONLY 

C  F 2  —  UPPER  CUTOFF  FREQUENCY,  PASS  AND  STOP  FILTERS  ONLY 

C  NT IMES  -  NUMBER  OF  DATA  POINTS 

C  VIN-  INPUT  DATA 

C  VOUT-  OUTPUT  DATA 

PI-4.0 DO »DAT4N*1.3DO) 

CALL  RDATA 
WC=2*DTAN*P1 *FC/FS» 

W 1 =P I *F 1 /FS 
W2-PI«F2/FS 

GO  TO  *10,20.301, ITYPE 
10  CALL  BUTTER 

tO  TO  40 
20  CALL  CHEB 

CO  TO  40 

30  CALL  CAUER 

40  CALL  FACTOR 

CALL  FPLOT 
JJ-ITRANSM 

GO  TO  (100, 50, 6 9, 70 1 »  JU 
50  CALL  HIGHP 

GO  TO  80 
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B— 1 .  MAIN  PROGRAM  (Cont'd) 


to  CALL  BANDP 

GO  TO  83 

73  CALL  BANOS 

83  CALL  F  PLOT 

133  lFINTlHES.EO.Olia  TO  110 
CALL  F  11  1HP 
CALL  P  L  T  DAT 
110  CONTINUE 

STOP 
END 


B-2.  SUBROUTINE  RDATA 


SUBROUTINE  RDATA 
REAL#8  PI 

COMMON  FC,  FS,  EPSI  ,  EPS?,  PI,  JTTPE,  N,  1  TRANS 
CDMMON/tfDATA/vmiOOO  1,  VOUT  11000 1  ,NT1MES 
COMMON/TRANS/  F  1  ,  r  2 

RE  AO  (S  ,1 0 1  HYPE,  N,  1  TRANS,  FC,FS,EPS1,  EPS? 

13  FORMAT (315,  4E1D.31 

F  2=  0 
F  1  =F  2 

IF  UTRANS.E3.0I  GO  TO  100 
READ«5,20»  F  1  ,F  2 
23  FORMAT  t?E10 .3) 

READ15,25)NTIMES 
25  FORMAT ! I  51 

IF  INTJHES.E0.01G0  IQ  50 
JJ= INT lMES«7}/8 
DO  100  K -1 , J J 
L  =  IK-1  1*8*1 
M*K  *8 

130  READf5,301IVINII 1,1 =L»M1 

33  F  ORMAT I8E10.31 

53  CONTINUE 

NR1 TEIb,45 ) 

55  FORMAT  I////  ,5X5tmYPE  ,7X  1HN  ,6XfcH  I  TRANS  .9X2HFC  ,  13X2HFS  ,  1 2X4HEP  S 1 , 
IUX4HEPS2,12X2HF1  ,1  3X2HF21 

VR1TEI6.401  ITTPE,  N,  JTRANS,  FC ,  FS,  EPSI.  EPS2 »  FI,  12 
43  FORMAT  1/ ,315X15) ,  b I5XE 1 0. 31  ,/ // / 1 

RETURN 
END 
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B-3.  SUBROUTINE  BUTTER 


SUBROUTINE  BUTTER 

R  E  A  L  *8  A  0  4  2  0  1  ,  41120),  BII20),  B2420),  PI, PH) 
C3MPLEX*16  P  I  20 1 t  Z 1 2  0  > 

COMMON  FC,  FS,  EPS1,  EPS2,  PI,  ITYPE ,  N,  1 TRANS 
COMMON/F 1LT/WC,  U 1 , U2  ,  AMP, AO  ,A 1 ,  A2  ,  61,  62,  P,  l 
REAL»8  A  MP , A  2 
A*N 

K=N-2* INTIN/2. > 

DO  20  lO.N 
B  =  I 

1FIK.EQ.0I  C 0  TO  5 

PHI -PI *4 (A«l .d)/A 

tO  TO  ID 

PHI *P1 #4 1 .0  /  I2.D*A) *IR-1 „0)/A«0.51 

Pill-  WC*OCMPIXIOCOSI PHI  I, DS IN (PH I  11 

CONTINUE 

AMP  =WC  **N 

RETURN 

END 


B-4.  SUBROUTINE  FACTOR 


SUBROUTINE  FACTOR 
RE  AC *8  AMP 

REAl*8  A0I2  ),  41  4  20),  B2I20I,  B2I20T,  PI,  RR , R 1  ,D ,  AK ,  6K,A2 

COMPLE  X*lb  P  4  20 ) ,  Z  420I,  R420),  X 

COMMON  FC,  FS,  EPS1,  EPS2,  PI,  1  TYPE ,  N ,  ITRANS 

COMMON/F ILT/KC,  U 1 , N?  ,AMP, AO ,A 1 , A2  ,  91,  62,  P,  Z 

A2=1.00d 

K*  I  NT (I N« 11/2.) 

DO  50  1=1, K 

X-DCNPLXC1.0D0,  0.0001 
00  20  J= 1 , N 
IF  4J.EQ.D  GO  TO  20 
X  =  X*  4P  1 1  I-P  I  J I) 

25  CONTINUE 

R ( 1 1 =  AMP/X 
IF4N.E0.1IG0  TO  33 
X  =  DCMPLX  41 .ODO.D.DDOI 
I F  4  I  TYPE .NE . 3  1&3  TO  50 
KK  =  2*IN/2) 

00  30  J=1,KK 
X  =  X*4PII  1-ZI  Jl) 
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B-4.  SUBROUTINE  FACTOR  (Cont'd) 


3D  CONTINUE 

R(I)=RII  )*X 
AMP-AMP*(KK+1-N) 

53  CONTINUE 

IF  Cl  TYPE .NE  .  3  )A1P=0 .000 
00  100  1  =1,  K 
RR=DREAL  (MI)) 

RI =0  JNAG  CR(  I  )  ) 

AK  =  -DREALJP  I  I  1M0.5DO 
BK  =  01  MAGI P (  I  )) *3 . 5  DO 
0=41. 003*  AK)**2*BK**2 
AOC I )= IRR*C 1  .0D3*AK  )  -RI*BK!/D 
Al( I ) =- I RR*  ll.33D-AK)*Rl =BK )/0 
BMI)=  -2.000*11  .  3DQ— AK**2— BK**2 1/0 
B24I)= (( 1.0DD-AIC»**2*BK**2»/D 
130  CONTINUE 
I =2*K-N 

IF!  l.EQ.OIGO  TO  113 
A0IK)=A3(K)/2 
B1 ( K )  =  B 1 ( K ) / 2 
A 1  {  K  )  =  0  • 

B  2  (  K  )  =  3  . 

110  RETURN 
END 


B-5.  SUBROUTINE  CHEB 


SUBROUTINE  CHEB 

COMMON/F  ILT/MC,  HI  ,  W2  , AMP, AO  ,A1  ,A2  ,  Bl,  B2,  P,  l 
COMMON  ft,  FS,  EPS1  ,  EPS2,  Pi,  I  TYPE ,  N ,  1  TRANS 
C0MPLEX*16  M  20 )  «  Z  420) 

RE  AL  *8  A  0 1 2  L  )  t  A1420),  B1I20I,  B2C20) 

REAL *8  CHV,SHV,Y,J,AMP,PI, A2 
EPS1 =SORT II .3/E PS  1* *2-1.01 
A2*l  .003 
A  =  N 

V  =  DL0G(1  .0  00/ EPSMD  SORT  II. 000*1  .QDO/E  PS1**2 1 )  /  A 
CHV  =  CDEXPIV  )  *DEXP»-VM/2  .0  00 
SHV =CDEXPtV>-OEXP I -V) 1/2.000 
JJ~  IN* 1 ) /2 
DO  50  1=1, N 

U=PI»I2.0DO*I1*N)-1 .ODOI/f 2.0D0*AI 

PI)  )  =NC*OCMPLXlDSIN (U)*SHV,  -DC 0S4U ) *CHV 1 
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B-5 .  SUBROUTINE  CHEB  (Cont'd) 


53  CQNT  INUE 

K- N-2«JJ 
AMP  =  1.000 
IffK.EQ.01GL  TO  93 
JJ=JJ-1 

IffN.EG.llGO  TO  S3 
DO  75  1=1, JJ 

AMP = AMP* CORE ALCP4 I  I  i • *2* D| MAC  I P I 11 1**2 1 
75  CONTINUE 

83  ANP  =  AMP*OR£AUPI JJ*11  > 

GO  TO  125 

93  DO  100  1=1, J J 

AMP=AMP*fDREALlPfll  I  **2*  DIMA6tPll))**2> 
130  CONTINUE 

AMP  =  AMP/DSOR  TI1.0D0-*EPS1  **2> 

125  RETURN 

ENO 


B-6.  SUBROUTINE  CAUER 


SUBROUTINE  CAUER 

COMMON/FILT/MC,  Ml , W2 , AMP, AO ,A 1 ,A2 ,  Bl,  B 2,  P,  l 
COMMON  FC,  FS,  EPS1,  EPS2,  PI,  1  TYPE ,  N,  1  TRANS 
REAL  *8  AMP,Q*A,S0,N  ,0  ,V 

REALMS  AO  1 20  > «  A1420),  B1I20),  82120),  PI,A2 
COMPLE  X*  16  2  1201  , PI20) 

S0=0.000 
M  =O.ODD 

G*C1 .003-08ie<EP52M*2>*M0.25» 

C=0. 500*441.000-01/ tl.ODO*0)> 

AMP  =  1 . 000 

G  =  0«2.0*Q**5* 15.0  *Q  **9*1 50 .0 *ft**l 3 
A  =  08LE IE  PSl  ) 

A =0 LOG  4  4  1 .000 /A* 1 .0  DO  1/4 1.000/ A- 1 .ODO 1 1 / 12 .0 *N1 
0=0.000 
OQ  10  1=1,100 
M  =  1  — 1 

SO=SO*  1-1  .ODO  1**M*0  **  IM*I1 *DSINH  1 4 I*M 1 *A1 
If (H.EQ.SOtGO  T3  ^0 
10  M  =  SO 

23  DO  30  1=1,100 

0=0* 1-1. ODO 1**1 •0**11 •1>*DC0SH|2.0*1*A» 

I F (O.EQ.M)GO  TO  93 

33  N  =  0 
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B-6.  SUBROUTINE  CAUER  (Cont'd) 


41  SO -DAB  SI  2. j*0®*!0 .2 5)* SO/ 11 .000*2.0*0 1) 

C  -N 

«=DSQRT|  II  .lDC*EPS2*50**2>*I1.0D0*S0**2/FPS2>) 
IFIN.E0.11G0  TQ  100 
J J  =  N/2 

U= I2*JJ* 1-N 1/2.3 
DO  100  1=1, JJ 
0  =  0. 0D0 

V  =  0 

DO  50  J*1»1C0«2 
M  =  J— 1 

0=1  -1.0)  ®*M*Q*MH*  J  I®DS1NI  CH*J1*PI*IFL0AT  II  1-U  1/0*0 
N=J*1 

U=|-1.01**J*Q**IM*J)®D5IN1  (H*J>*Pl*fFLOAT  II  1-Ul/C1*Q 
IF  fO.EQ.VIGO  TO  63 

53  V  =  0 

63  A=0 .ODO 

V  =  A 

DO  70  J=  1  (  1  0 0  >2 

M  =  J 

V  =  V*l-1.01**N*Q®®IH®f1  l®DC0St2.0®«®Pl®IFL0AT»II--U)/C» 
K  =  J*I 

V=V*I~1.0»®*M*Q**IH®M 1*DC0S»2.0*M*P1 *1  FLOAT! 1 1  — U I/O 
1FIV.EQ.A1GO  TO  83 
73  A  =  V 

80  0=2. 0*0* *10 .251® 0/1 1.0DO*2.0*V I 

V  =  I 1.0  03-EPS2*0®OT®4J .000— 0*0/ EP 52  I 
V=DSQRTIV) 

Zll  l=DCNPLXt 0.030,1  .0  00/0 1  *i<C 
Z •  N ■*  1—  1  1  =DCGNJGt  Zl  1  11 
A=  1  .ODO*  I  SO  *  0 1**  2 
P«11=DC!1PLX*-50®V/A  ,D*W/A1*MC 
PCN*l-ll*DCONJG{Pti >1 

AHP=AMP*0*0*  f  CS3*Vl *«2*I0®K1 «®21/A**2 
130  CONTINUE 

K=N-2*«N/21 
I  FIK  .EQ.OIGG  TO  150 
PIIN*11/ 21  =  DCMPL  X I -SO  *MC ,0 .ODO 1 
AHP  =AMP®S0*NC 
GO  TO  200 

150  AMP=AMP®EPS1 
230  CONTINUE 
RETURN 
END 
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B-7.  SUBROUTINE  HIGHP 
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SUBROJT I NE  HIGHP 

RE Al*B  AOt  2  »*  H 12  0)  •  B1I201,  B2I20),  PI 
REAL*8  A  MP , A  2 

RE  A  L  *8  AP.Ai  D,A11,B11  ,B22,D 

COMMON  FC,  FS,  cPSl  ,  EPS2,  PI,  1  TYPE ,  N,  1TRANS 
COMMON/F ILT/WC,  W  1  , M2 , AMP , AO , A 1 , A2 ,  Bl,  B 2,  P,  l 
M3=P1*FC/FS 
A2S-1.0D0 

AP*-casiui«W3i/:as(ui  -M3» 

DO  100  K=1,N 

D  =  1 .0D0-81IK  >*AP*B2  IK  J»AP*AP 
AQU-IA0IKI-A1 (K1 * AP >M1  .ODO-AP  l/D 
All  *IA0IK)*AP-U  IK)  J« 11. ODO-AP 1/D 
B11MAP*AP-B1  IK  I  *41  ,ODO«AP  *AP  I«2  >0D0*B2  CK 1  *API/D 
B22=  (AP*AP-AP«=BHR)  +B2IKH/D 
AOlK)=AOO 
A1 f K)  =  A1 1 
Bl ( K I s  B 1 1 
B2CK I-B22 
130  CONTINUE 
RETURN 
END 


B-8.  SUBROUTINE  FPLOT 


SUBROUTINE  FPLOT 

COMMON  FC,  FS,  EPS1,  EP52,  PI,  1TYPE,  H,  ITRANS 
COMMON/F ILT/MC,  Ml ,M?  ,AMP,AO,AI ,A2,  Bl,  B2,  P,  I 
REA L*8  AMP 

REALMS  A0C201,  A14201  ,  B1I20I,  B2I20),  PI.A2 
DIMENSION  F 12511 ,HE  JMTI251 1, PHI  1 251 1 
DIMENSION  X  LABI  2) ,Y LAB1C2I , YLAB2 12 1 ,HE AD«1 0 1 .SUBHI21 
DATA  YLAB2I1  ),YLAB2I2I/8H  PHASE  ID  ,BHEG 1  / 

DATA  YLAB1I1 ) ,VLAB1 12 1/8HAMPL 1TUD ,BHE  / 

DATA  XLABI1I.XIABI2 )/ 8HF RE QUE NC , 8HV  fF/FSl/ 

MR  I TEI6.800  1 
JJ= IN* 11/2 
DQ  805  I*1,JJ 

8  35  MR  I TE 1 6, 801 1A0II I, All  1), Bill  I, 82(11 

MR1TEI6, 802)4 MP,A2 

830  FOR  MAT llbX2HA",32X2HAl ,32X2HB1 ,32X2H82 ,/// 1 

831  FORMAT  II P4 ( 1 XQ32 >25 ) ) 
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B-8.  SUBROUTINE  FPLOT  (Cont'd) 


0 02  FORMAT  I4X6HAMP  =  ,1PQ3?.?5,10X5HA?  =  ,0PF4.0I 
JJ  =  IN* 11/? 

READ«S,10)  i HEAD  II), 1=1, 10) 

13  FORMAT  1 1  0A8 ) 

DO  100  1*1,250 
A=  1-1 

Fill  -A/249  .<• 

U  =  2.0*Pl *F(  1  1 
C=CQS(U) 

S-S 1 N t  W 1 
PHI  411=0  . 

HE JWTI 1) =0 . 

DO  50  K  =  1 ,  J  J 
CN= AOIK) *A1  IK)*: 

DM*— A 1  IK )*S 

C  0  =  1 -0  *B 1 1 K ) *C*S2 IK !• 42.0*C*C-1 .0 1 
DD=BltK)*S*2.0*52IK )*C*S 
DD= -OD 

EE=CO*CD*DD*DD 
FF=4CN*CD*DN*DD)/EE 
E  E  = I DN*C  D-DC  *CN 1 / EE 
HE JMTI1 1 =HE  JMTII l«FF 
PHI  II) =PH1«  1)*EE 
50  CONTINUE 

FF  =  |1.0D0*A2’»CI*HE  JWT 1 1 > «S»PH1 1 1 ) * A2  *AMP 
EE  =  H.0D0*A2*CI*PHI 11  )-S*HE JNT I I) *A2 
HEJWTtl) *SQRT IEE  *EE  *F  F  *F  F 1 
PHI  II I  * ATAN2  IEE  , F F  1 
PH  1 1 1 ) *PHI 1  I  IHJO.O/P  I 
100  CONTINUE 

CALL  DRAN1DI1 ,4,4 ,20,0,250 ,2.,0.,XLAB, 

1  YLAB1,HEAD, SUBH.F ,HE JWT ) 

Fll  1*0.0 

CALL  ORA N10 1 1  A ,20,0,250 ,2.0 ,0.,XLAB, 
1YLAB2, HEAD, SUBH,F, PHI  I 
RETURN 
END 
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SUBROUTINE  BIND? 

CQHHON/F 11T/WC.  U l , M2  ,AMP,A0 ,A1 ,A2 *  81,  B2V  P,  l 
COMMON  FC,  FS,  EPS1,  EPS 2.  PI,  ITTPE,  N,  ITRANS 
REAL  *8  AOf  20 1  ,  A1I20),  BK20I,  B2I20),  PI,A?*AMP 
COMPLEX*  16  2  120)  ,P(20),PP,0,R,S,T,U,V 
RE AL*8  C  ,D,AL  »E,F 
A2=Q .003 

AL-OCOS(DBLE(WI*W2I l/DCO  SID8LEIW2-RII) 

JJ=N/2 

IFIN.EQ.1IGO  TO  S3 

DO  50  1=1, JJ 

AMP  =AMP  * AO  I  I  1 

C=B2fl)-81ll)**2/4.0D0 

lFIC.LE.O.OOOlGO  TO  40 

PP  =  DCMPLXI-B 1 (I ) / 2 .000 ,DSQRT 1C ) 1 

Q=fA0fl)*PP*PP*CA3II)*Al II ))«PP«OCMPLX (All 1 > v0. 000)1 
1/DCMPLX|O.ODO,2.030*D1MA6IPP1) 

V*DCMPLXC1. 000*3. 3D0I 
S  =  AL*IV«PPI/2  -ODO 
R=S*CDSaRTI5*S-PP| 

S=S-CDSaRT«S*S-PPI 
T*0* IAL* R-V I / IR-S I 
U=Q*tAL*S— V1/IS-R1 
C=DREALIR) 

D  =  D  IMAGIR) 

E  =DREAL IT) 

F  =  D IMAGI T| 

F=-2.0*IC*E*0*F) 

o=c*c*o«o 

C=-2.0*C 
E=2.0*E 
AMP  =AMP*F/D 
AO  1 1  )  =  — F/0 
A 1 ( I  )  =  E-F*C/D 
Bl 1 1  )  =  C 
B2(  I  )  =  D 
K= I N*1 1/2*1 
C=DREAL(SI 
0=0 1MAGIS) 

E  =DREAL(U) 

F=D IMAGIUI 

F=-2.0*IC*E«0*f  I 

D*C*C«D*D 

C  =  -2.0*C 

E*2.0*E 

AMP  =AMP«F/0 

A0«KI*-F/0 

A]  IK)=E-F*C/0 

B1IK)*C 

62<K)*D 

GO  TO  50  59 
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B-9.  SUBROUTINE  BAN DP  (Cont'd) 


43  C -D  SORT  I  -C  I 

D*-Blll>/2.  jDO*C 
C=-B1I1»/2.0D0-C 

E  =  IA0I1)*C*C*  IA3I1 )«A1I1  1)*C«R1I1)>/«C-0I 
F= I AOI 1) •D*D« IR3 1 i )  «A  1 II ) } *D«A 1 II) 1/ID-CI 
AMP=AMP-E/C~F/D 
AO  I  1 1-E/ C 
All  I 1=-AL*£/C 
B 1 1  1 )  =  -AL»t 1 .0D3*C1 
B2|I1»C 
K*(N*1 1/2*1 
AOIK  l-F/D 
A1IK)*-AL*F/D 
B1«KJ=-AL*I1.0D3*D1 
B2 1 K  )SD 
53  CONTINUE 
K*N-2»JJ 

I F IK .EQ.01GQ  TO  130 
AMP  =AMP« AO  INI/BIIN) 

D-Bl  IN  1 

C*A0INI*I1.  CD  0-1 .3/01 
K*IN«l)/2 
AO  IK  MC 
A1 I K )*-AL*C 
B1|K)*AL •ID-1. 0301 
B2IKH-D 
130  N=2*N 

RETURN 
END 


B-10.  SUBROUTINE  BANDS 


SUBROUTINE  BANDS 

COMHON/F ILT/tfC,  N 1 »N2  ,AMP, AO ,A 1 , A2 ,  Bl,  B2.  P,  Z 
COMMON  FC,  FS.  EPS1,  EPS 2»  PI,  ITYPE,  N,  ITRANS 
REAL *8  AO  120 )  ,  A1I201  ,  81(201,  621201,  P!,A2,AMP 
COMPLEX*  16  ZI20),PI201,PP,QtR,S,T  ,U,V 
REALMS  C,D,AL,E,F,X,C1,C2 
A2*0 .000 

AL  *DC05 I DBLE ( W1 *W  2 ) 1/DCOSIDBLE IN2-N1 1 1 
JJxN/2 

X*OTANIOBLE IN2-A11I*DTANIP1*FC/FS» 

C1*2.0*AL/IX«1.3031 

C2r 1 1 .ODO-X 1 / IX* 1 .0  DO  1 
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B-10.  SUBROUTINE  BANDS  (Cont’d) 


IFIN.EQ.DGC  TO  5 ) 

DO  50 

v=dcmplxh.odo.3.ddoi 

AHP  =  AHP*  AO (  I  ) 

C=B2tIl-Bl(I) *•?/%- 00 0 
IMC.LE.0. 000)03  10  AO 
PP=DCMPLXt-Bl  I1)/?.0D0*D SORT  1C  11 

Q-(AOIJ)  *PP*PP*IA0t  ))  «A1  I!  »»«PP*DC»PLXIAlin  rO.ODOIT 
1/DCMPl XI 0.000,2. ODD *D1MAG(PP»» 

Q=Q*C2/IV-PP*C2) 

AMP=ANP«2.O*0AEALI  0) 

S  =  C1»IV-PP)/U¥-C?*PP  )*2 .0  > 

T=tC2*V-PPi/{V-:2*PP) 

R=S*COSQRT(S*S— T  1 
S  =  S-C0SQRTI5«=S-r  J 
T=Q*«R*R-C1/C2*A*V/C2I/CR-S) 
U=Q*(S*S-tl/C2*S*Y/C21/IS-*} 

C=DREAl ( R } 

0=0 IMAG ( R ) 

E  =  DREAL 17) 
f-DIMAGIT) 

F  =  -2.0*(OF  > 

0=C*C«D*0 
C=-2.0*C 
E  =  2.0*E 
AHP  =AHP*F/0 
AOI 1»*-F/D 
AUI  )=E-F«C/D 
Bltl>=C 
B2d  )SD 
A  =  I  N-*l  i/2*  I 
ODREAUS) 

0=0 IHAG(S) 

E  =DREAL  tU) 

F*01KAG(U» 

F  =  -2.0*IC*E-*D'>F1 

D=C*C«D«D 

C  =  -2  .0*C 

E*2 .0*E 

AHP=AAP*F/D 

A0IK)=-F/0 

AllKI=E-F*C/0 

#1  m  >=c 
62(10  =  0 
GO  TO  50 


\ 
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B-10.  SUBROUTINE  BANDS  (Cont'd) 


43  C  =D SORT  l-C ) 

D*-BllJ)/2.0Dd*c 
C=-B1I  11/2.0D0-C 

EMA0in*C«C»U3(lHAnill*CfUm!/iC-Dl 
FMAOt  I)  *0«D*  «AD  (  I  ) *A 1 II) ) *D*A1 ( I > )/ t D-C ) 
AHP=AMP*E/IC2-C)*F/  1C  2-0  » 

AOI I )=E«IC2/C1.3D3«C2 1-1 .0/IC2-C)) 

A 1 1 1 1 s -E  *C  1  *  IC«C2 )  /II  .0D0«C2)**2 
Bit  1)=C1*IC-  1.030)  /II  -ODO*C  2  ) 

B2( I »  =  €  C  2-C >/(1.0D0«C2) 

K=I«N 

B2 IK )s (C  2-D )/ll*000«C2) 

B1 I K)*C1*ID-1  .0001/11  .0D0«C2) 
A0lKJ=F*IC2/|I.3D3«C2)-l  -0/IC2-0H 
A 1  lK)=-F#CI*fD«C2>  /  II  .0D0*C2  1**2 
SO  CONTINUE 
K  =N~2* JJ 

IF(K.EQ.O)  GO  13  100 
K=lN«U/2 
AMP=AMP«A0|K1 
D=B1  IK) 

C*A0CKI/I1.0D0«B1IK ) *C2> *1 l .ODO-D ) 
B2lKl=IC2*D  )  /  II .  3  30  *C 2*0  I 
B1  IK)*-C1*ID*1.DD31/C1.0D0«C2*D1 
A HP  =AMP«C/B2 IK) 

A0IK)=C*C2~C/B2IK) 

A 1 (K ) *— C 1*C -  C*il IK 1/B2IK) 

100  N*2«N 
RETURN 
END 


B- 11.  SUBROUTINE  FILIMP 


SUBROUTINE  FILHP 

COMMON  FC,  FS,  EPS1  ,  EP52,  PI,  I  TYPE ,  N,  1  TRANS 
REAL*8  AMP 

REALMS  AO  1 20 )  *  UI20),  B1I20),  B2I20),  PI.A2 

COMMON/tfOATA/  VI N 1 1 00 0 1, VOUT 11000), NT IMES 

COMHON/F 1LT/NC,  M 1 , M2 ,AHP, AO ,A1 , A2 ,  Bl,  B2,  P,  Z 

RE A L*B  YM2.2C),  YSUM1,  YSUM2,  YTEM 

TSUM1*0.0 

YSUM2s0 . 0 

JJ-(N«l)/2 

00  10  K : 1 • U J 
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B-11.  SUBROUTINE  FILIMP  (Cont'd) 


YM2,K»  =  A0(K)*VINI1  ) 

VKf  1  rK)  =  A0{K>*VIN<2H<Al<K)-BllK)»A0tK>l«VlNCl> 
YSUH2sYSUM2+VK(2*K) 

YSUH1  =  YSUH1-»YMI  ,<> 

13  CONTINUE 

VOUT  II  )  =  Y5UM2 
V3UT 12 1  *  YSUH 1 
DO  50  J*  3»NT I  HE  S 
YSUH2--YSUH1 

rsuNi*o.o 

DO  40  K  *  1 »  JJ 

YTEN=AOIKl*VlNtll*A 1 IK )*V1N  t J-l 1 
1-lKKinMl  ,K  1-02 IX  1*YKI2,R) 

YKI2  »K 1*  YK 1 1 »R) 

YK  1 1  »K  1  *  Y1  EM 
YSUM1*YSUM1*YTE1 
43  CONTINUE 

VOUT  I J)  -Y5UN  1-»VSUN2  *A2  ♦AMP«VlNf J) 

53  CONTINUE 

RETURN 
END 


B-12.  SUBROUTINE  PLTDAT 


SUBROUTINE  EL  TO AT 

CONMON/YDAT A/ VISE  V 300  1  .YOU! 1 1000  1 , NT  I  ME S 
COMMON  FC.FS 

DIMENSION  NE«DIIQ>,XIA8I2),»LABI2|.TC10001 
DATA  XlABtl>'XUBI?)/BNTtME  ISE  » BHCONDS 1  / 

DATA  SUB1.SU82/SNV1N  .BHYOUT  ✓ 

DATA  YLABI1) ,  YUB 1 2  1/8HANH1TUD  ,BHE  / 

DO  20  1*1»NT 1MES 
23  T( 1 1 *( 1- 1 l/f  S 

READ (5 *10) « HEAD* 1), 1- 1,101 
13  FORMAT  1 1 OAB  1 

CALL  0RAN10I1  ,4,4 ,20, 2tN TIME S, 2.,0 ..XL AB, YCAB , HE AD ,SUB 1 ,T , V1N 1 
READ (5. 101 (HE ADI  11,1-1,101 

CALL  0RAN10I1  ,4, 4 ,20,2 ,NT1ME 5,2., 0.,XLAB*YLAB, HE  AD, SUB 2.T, VOUT’ 

RETURN 

END 
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